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Abstract 



^SJ ' Motivated by the negative thermal expansion observed for sihcon between 

' 20 K and 120 K, we present first an ab initio study of the volume dependence of 

' interatomic force constants, phonon frequencies of TA(X) and TA(L) modes, 

o ■ 

lO ' and of the associated mode Griineisen parameters. The influence of successive 

nearest neighbors shells is analysed. Analytical formulas, taking into account 
interactions up to second nearest neighbors, are developped for phonon fre- 

> ■ 

S ' quencies of TA{X) and TA(L) modes and the corresponding mode Griineisen 

> 

' parameters. We also analyze the volume and pressure dependence of various 

■ thermodynamic properties (specific heat, bulk modulus, thermal expansion), 

and point out the effect of the negative mode Griineisen parameters of the 
acoustic branches on these properties. Finally, we present the evolution of 
the mean square atomic displacement and of the atomic temperature factor 
with the temperature for different volumes, for which the anomalous effects 
are even greater. 

PACS numbers: 65.70.+y; 65.50.+m; 65.40.-f; 63.70.+h; 63.20.Dj 
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Typeset using REVT^ 



INTRODUCTION 



In the past few years, theoretical and algorithmic advances have made it possible to de- 
termine the thermodynamical properties of solids (such as specific heat or thermal-expansion 
coefficient) from first principles calculations 

Silicon is a choice system for testing these methods, since accurate measurements on 
high-purity samples exist over a wide range of temperature. Moreover, it presents a nega- 
tive thermal-expansion coefficient between 20 K and 120 K which is of fundamental interest. 
This unusual thermal-expansion behavior can be attributed to the negative Griineisen pa- 
rameters (i.e., phonon frequency increases as crystal volume increases) of the transverse 
acoustic (TA) phonons near the Brillouin-zone boundary. It is interesting to note that this 
anomalous behaviour also affects other properties of silicon, such as the mean square atomic 
displacement (as shown later in this paper). 

This anomalous negative thermal expansion has been analyzed in previous studies, most 
of them relying on the quasiharmonic approximation 0. In this approach, the volume 
dependence of phonon frequencies must be determined. Kagaya, Shuoji, and Soma 0] used 
a perturbation treatment with a model pseudopotential to calculate the specific heat and 
the thermal-expansion coefficient. Biernacki and Scheffier employed a Keating model 
with two parameters, which were extracted from local-density pseudopotential calculations, 
to compute the thermal-expansion coefficient. Fleszar and Gonze |^ performed a model-free 
computation of the thermal expansion, using a linear-response technique within DFT . Xu 
et al. 1^ used a tight-binding model to calculate the thermal-expansion coefficient and mode 
Griineisen parameters; they also showed by using a simple model that the negative values of 
the mode Griineisen parameter could be attributed to a larger contribution from the central 
part of forces than the angular part of forces. More recently, Wei, Li, and Chou extracted 
interatomic force constants from ah initio calculations of planar forces and calculated 
the specific heat, the overall Griineisen parameter, and the thermal-expansion coefficient. 

In this paper, we use a variational approach to density-functional perturbation theory 
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to calculate volume-dependent dynamical properties of silicon (see Sec. |). We present first 
an ab initio study of the volume dependence of interatomic force constants up to twenty- 
fifth nearest neighbors. Wei and Chou had also presented such a calculation but only for 
one volume and up to eighth nearest neighbors |TO|JTl[] . Phonon frequencies of TA(X) and 
TA(L) modes, and of the associated mode Griineisen parameters are also calculated for 
different volumes. The influence of successive nearest neighbors shells is analysed. We 
confirm that the contribution of atoms connected by zig-zag chain along [1 1 0] direction 
are the most important, as suggested by Mazur and Pollmann |jl2|. But we show that the 
contributions of flfth, sixth, and seventh atoms along the chain (respectively thirteenth, 
eightenth and twenty-fifth nearest neighbors) are not negligible contrarily to what had been 
speculated by Wei and Chou [|10[- Analytical formulas, taking into account interactions up 
to second nearest neighbors, are developped for phonon frequencies of TA{X) and TA(L) 
modes and the corresponding mode Griineisen parameters. In Sec. |T[ we analyze the volume 
and pressure dependence of various thermodynamic properties (specific heat, bulk modulus, 
thermal expansion). The effect of zero-point motion on static equilibrium properties (lattice 
constant and bulk modulus) also analyzed. We point out the effect of the negative mode 
Griineisen parameters of the acoustic branches on these properties. In Sec. we present 
the evolution of the mean square atomic displacement and of the atomic temperature factor 
with the temperature for different volumes. Anomalous effects present at all temperature 
are investigated, using a band-by-band decomposition. Finally, we present our conclusions 
in the last section. 

I. DYNAMICAL PROPERTIES 

A. Interatomic force constants 

To obtain accurate interatomic force constants (IFC's), we first calculate from first prin- 
ciples the dynamical matrices for 10 different wavevectors in the irreducible Brillouin zone 
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(IBZ) (the special points mentioned in Ref. [T^) using a variational approach to density- 
functional perturbation theory Then, we perform a Fourier Transform of these dynamical 
matrices in order to get the IPC's, following Ref. This sampling of dynamical matri- 

ces allows to obtain IFC's in a real space box containing 512 atoms. This is more accurate 
than what had been done in Ref. where the IFC's had only been calculated for 2 different 
wavectors in the IBZ which corresponds to a real space box of 64 atoms. It is also better 
than the work of Wei and Chou [0,0] who had included interatomic interactions up to 
eighth nearest neighbors, which corresponds to 99 atoms. 

Our calculations, performed within the local density approximation (LDA) [T^, use a pre- 
conditioned conjugate gradient algorithm |1T7|JT8[| . We use a rational polynomial parametriza- 



tion of the exchange-correlation energy functional |jT9[, which is based on the Ceper ley- Alder 



gas data The electronic wavefunctions are sampled on a mesh of 10 special k points 
in the IBZ and expanded in terms of a set of plane- waves whose kinetic energy is limited 
to 10 Hartree. The "all-electron" potentials are replaced by an ab initio, separable, norm- 
conserving pseudopotential built following the scheme proposed in Ref. PT|. The calculated 



equilibrium lattice constant is 10.18 Bohr (the influence of the zero-point motion on this 
value will be discussed in Sec. |Tp, whereas the experimental one is 10.26 Bohr. 

The calculation of IFC's is performed for three different volumes, corresponding to lattice 
constants a of 10.00, 10.18, and 10.26 Bohr respectively. The IFC's for these different lattice 
constants are listed in Table |, in which coordinates are in units of a/4 and notations for 



the force constants follow Ref. p3 



As already noticed in Refs. IT0IJT4 , IFC's for atoms connected by zig-zag chain along 
[1 1 0] direction (or the directions related to it by symmetry) are the most important. The 
magnitude of IFC's along these peculiar directions decays quite slowly. The contribution of 
the first nearest neighbors (ai and jSi) is about 4 to 5 times smaller than on-site contribution 
(ao). Going up to second nearest neighbours reduces the IFC's by another factor of 5 (A2) 
or 10 (/i2 and 1^2). The contribution A25 of the seventh atoms along the chain (coordinates 
(7,7,1) in units of a/4), which are already the twenty-fifth nearest neighbors, is still of the 



order of 1 percent of the contribution of the first nearest neighbors, and of 5 percent of that 
of the second ones. Comparatively, the biggest contribution of the seventh nearest neighbors 
(these do not belong to the chain, but are twice closer than the seventh atoms along the 
chain) is about 2 times smaller. 

In order to vizualize these IPC's, we compute the forces induced along the zig-zag chain 
by the displacement of a generic atom of this chain in three perpendicular directions (along 
[110] — parallel to the chain direction, along [0 1] — perpendicular to the chain direction, 
but in the plane of the chain, and along [110] — perpendicular to the chain direction, but 
out of the plane of the chain), all the other atoms being kept fixed. These forces can easily 
be obtained from IFC's by: 

F^o. = -$a/3(«:; k') {R - Ro)^,p (1) 

where R — Rq is the displacement of the generic atom, and the IFC ^^/^(/t; k') is the force 
exerted on ion k in the direction a due to the displacement of ion k' in the direction /3. 

The resulting forces have been reproduced in Fig. |^. The directions of the forces induced 
by the displacement of the generic atom in the [1 10], and [0 1] directions show that 
the bonds along the chain will tend to bend in order to keep constant the bond angles : 
at sufficient distance, the forces are nearly perpendicular to the last connecting bond. The 
decay of the magnitude of these forces along the chain is similar to that of the IFC's. At 
the contrary, the forces induced by the displacement in the [1 1 0] direction are much more 
smaller. Indeed, they have the generic form (/x — z/, z/ — yU, 0) where /z and v are of the same 
order of magnitude (see Table |l|) . It is thus more difficult to interpret them. 

B. TA(X) and TA(L) phonon frequencies and associated mode Griineisen parameters 

In previous studies, the negative mode Griineisen parameters associated with the modes 
of the acoustic branches, near the zone boundaries, has been shown to drive the negative 
thermal expansion. We now focus on the analysis of the high symmetry TA(X) and TA(L) 
modes because analytical results can be obtained for them. 
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The phonon frequencies can be calculated by solving the dynamical equation: 

"^api^; k')u„{k'P) = M^uolu„{Ka) (2) 

where Ua{i^a) is the displacement in the direction a of ion k for the normal mode a, is 
the mass of ion and uja- is the frequency of the normal mode a. 

In order to analyse the influence of the successive nearest neighbors shells for TA(X) 
and TA(L) modes, we artificially limit the interactions taken into account in the IFC matrix 
k') to the atoms whose distance is less than a cut-off radius. Then, we increase the 
cut-off radius shell by shell until the sphere contains all the atoms included in the real space 
box defined by our sampling of the Brillouin zone. 

In truncating the interaction to different shells, we break the acoustic sum rule (which 
states that moving all atoms as whole should not create any force) unless it is reimposed 
by artificially changing the on-site IFC a^. We analyzed both the results obtained with and 
without reimposing the acoustic sum rule and found no significative difference regarding the 
convergence with respect to the number of nearest neighbors shells included. Thus, in the 
following, we will only reproduce the results obtained without reimposing the acoustic sum 
rule. 

The results are reproduced in Tables || and |TTI[ The well-converged calculated frequencies 
of TA(X) and TA(L) modes are respectively 140.46 cm~^ and 108.626 cm^^ at the calculated 
equilibrium lattice constant, 147.37 cm~^ and 112.86 cm~^ at the experimental one. These 
values are in very good agreement with experiment |2^: 149.77 cm~^ and 114.41 cm~^. 



The associated mode Griineisen parameters are calculated by: 

where a is the index of the mode, and V is the molar volume (A^a x volume of the unit 
cell), linked to the lattice constant by: 

V = AT^^. (4) 



6 



The results, reproduced in Tables H and |T|, are very sensitive to the volume variations. The 
calculated mode Griineisen parameters of TA(X) and TA(L) modes are respectively -2.295 
and -1.814 at the calculated equilibrium lattice constant. These values are not in very good 
agreement with experiment ||2^: -1.4 and -1.3. However, if we consider the values -1.782 
and -1.451 obtained at the experimental equilibrium lattice constant, the agreement is much 
better. 

Considering for reference the values obtained by taking into account all the atoms in- 
cluded in the real space box, we see that the contribution of atoms connected by zig-zag chain 
along [1 1 0] direction are the most important for the convergence of the phonon frequencies 
of TA(X) and TA(L) modes, and of the associated mode Griineisen parameters. Beyond the 
fourth atom along the [1 1 0] chain (eighth nearest neighbors), the contribution of succes- 
sive nearest neighbors not belonging to the chains is not significant. But the contributions 
of fifth, sixth, and seventh atoms along the chain (respectively thirteenth, eightenth and 
twenty-fifth nearest neighbors) are not negligible, contrarily to what had been speculated by 
Wei and Chou |T^. Including still further atoms along the chains would necessitate to use 
another sampling of the dynamical matrix wavevectors in order to get a bigger real space 
box. Tables y and ^ also show that, once interactions with the first nearest neighbors are 
included, the anomalous behaviour of TA modes at X and L points (negative Griineisen 
parameters) are reproduced. This reinforce the pertinence of the model proposed in Ref . . 

For TA(X) and TA(L) modes, the eigenmodes U(j{K,a) are known for any ion from sym- 
metry considerations (see Fig. [l|). It is thus even possible to obtain an analytical expression 
of phonon frequencies, from the knowledge of the IFC's. 

If we only consider on-site interaction (just atom k) in the sum of the left-hand side of 



for both TA(X) and TA(L) modes, where m is the mass of silicon ion. The associated mode 



Eq. (D, we get: 




(5) 
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Griineisen parameters can easily be obtained by inserting this result in Eq. (Q): 

7 = 7(ao) (6) 

where 7(00) is the "force Griineisen parameter", that we define as: 
Icilnao 

^^"°) = -2dhn7- 

The force Griineisen parameters have been listed in Table H for for all IFC's whose magnitude 
is higher than 10"*^ Hartree/Bohr^. 

If we consider interactions up to first nearest neighbors atoms, we get: 



ujta{x) = y (8) 

for TA(X) mode, and 



^TAiL) = \J (9) 

for TA(L) mode. In fact, whatever the number of nearest neighbors taken into account, 
the phonon frequencies will always be written as the square root of the ratio of a linear 
combination of the IFC's and the mass. The associated mode Griineisen parameters can 
be written as weighted sum of the force Griineisen parameters corresponding to the IFC's 
under the root in Eqs. (||) and (|^) . For example, for the TA(X) mode, inserting Eq. (|^) in 
Eq. (D, we get: 



_ V d^ap + 4:(3i _ 1 V dao + 4/?i 

^^^^^^ ~ Vao + dV ~ 2ao + iPi dV ' ^ ^ 

This can be rewritten as follows: 

^^^^^^ 2\aQ + APij a^dV 2 I ao + 4/?i j p^dV ^ ' 



The definition of force Griineisen parameters, Eq. (|^, leads to 



iTAix) = (^\^] 7(«o) + (^^] 7(A) (12) 
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for TA(X) mode, and 

1TA{L) = ( , ^ ^r, ] 7(Qo) + ( , , ] 7(ai) + 

\ao + 2ai + 2(3i I \ao + 2ai + 2Pi I 



2(3 



1 



ao + 2ai + 2pi 



7(/3i) (13) 



for TA(L) mode. 

Interestingly the same argument holds, whatever the number of nearest neighbors taken 
into account. And thus, the same kind of formulas will be obtained for phonon frequencies 
and associated mode Griineisen parameters. 

Indeed, if we consider interactions up to second nearest neighbors atoms, we get: 



/ ao + 4/?i -4A2 

^TAix) = y (14) 

for TA(X) mode, and 



/ao + 2ai + 2/?i -4z/2 

^TA{L) = \j (15) 

for TA(L) mode. The associated mode Griineisen parameters can be written as weighted 
sum of the force Griineisen parameters corresponding to the IFC's under the root in Eqs. ([T^ 
and dll): 

-4A2 



V ao + 4/^1 - 4A2 / V "0 + 4pi - 4A2 / 



\ao + A/3i -4A2^ 
for TA(X) mode, and 



7(A2) (16) 



'^'^'^ = (ao + 2a,+V-4.2) ^^"^^ + (co + 2a,?2A - 4.2) ^^"^^ + 

^ ^ 7(A) + ( , . ~^"L . ] 7(^.) (17) 



\ao + 2ai + 2Pi - 4:U2 J V"o + 2ai + 2/5i - 4z/2^ 

for TA(L) mode. 

It should be noted that the weights in Eqs. (|T^ to (^) are not necessarily included 
in the interval [0,1] and can thus be negative. This explains why, though almost all force 
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Griineisen parameters are positive (see Table |), a negative mode Griineisen parameter can 
be obtained. 

For example, in Eq. (p!2|), we get for a=10.18 Bohr: 

-fTAix) = 3.08 X 1.1034 + (-2.08) x 1.9817 = -0.72 , (18) 

which is just what is obtained in Table 0. If all force Griineisen parameters were equal to 
1, we would also get 1 for the mode Griineisen parameter. This shows that the origin of the 
negative Griineisen parameter is the rather important difference between 7(/3i) and 7(0^0) • 

II. THERMODYNAMIC PROPERTIES 

Using the calculated phonon frequencies, the temperature-dependent phonon contribu- 
tion AF to the Helmholtz free energy F{V,T) is calculated as described in Ref. 0. This 
contribution is added to the energy of the static lattice E(y), calculated previously to 
get F{V, T) for a set of three volumes and for various temperatures. This function is then 
interpolated as a function of by a second order fit. The equilibrium volumes (or lattice 
constants) at various temperatures are determined by minimizing F{y^ T) as a function of 
V , as shown in Fig. ^. From F{y, T) any other thermodynamical property can be accessed. 

We first analyze briefly the specific effect of the zero-point motion. We find that the 
zero-point contribution to Helmholtz free energy F is AFo=12 J/mol, which is about 2.5% 
of the cohesive energy (460 kJ/mol p6[). This zero-point contribution causes the lattice 
constant to be shifted from 10.1894 Bohr to 10.1974 Bohr, which is a change smaller than 
0.1 %, and the bulk modulus Bt from 1.0387 Mbar to 1.0292 Mbar, which is a change of 1 %. 
Note however that the change in the bulk modulus includes two effects: first, it is linked to 
the second derivative of FiV) that includes the zero-point motion contribution AFq; and 
second, it is calculated for a different volume due to the shift of the lattice constant. In 
Fig. ^, we present the temperature dependence of Bt for different pressures. 

The entropy can be calculated following Ref. §. We get 5(298.15 K)=19.3 J K^^ mol'S 
to be compared to the experimental value of 18.81 J mol~^ We have also calculated 
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the variation of enthalpy H = F + TS - PV between and 298.15 K. We get 3.285 J mol"^ 
whereas the experimental value is 3.217 J mol~^ [^. We present in Fig. |^ the molar 
constant-pressure specific heat Cp^m for various pressures. This result agrees quite well with 
experimental work [^,^]. It is interesting to note that, for temperature higher than 85 K, 
Cpm is higher for low values of the pressure. Whereas, for temperature lower than 85 K, it 
is just the contrary. This means that there exists a temperature around 85 K for which Cp^m 
is independent of pressure. In order to understand this observation, we use the following 
relation (see Ref. ||30i): 

(19) 



'dS\ _ _ /dv' 



to write: 



This shows that in order to have Cp^ independent of pressure, one must have a decreasing 
ap. This is precisely the case of silicon between 20 K and 100 K, where the thermal-expansion 
coefficient gets more and more negative. 

This thermal-expansion coefficient is well reproduced by our calculations (see Fig. 
It is worth noting that increasing the pressure reinforces the anomalous negative behaviour 
of this property. This is confirmed by the strong negative value of the overall Griineisen 
parameter at high pressure (see Fig. This anomalous behaviour has also consequences 
on other properties. We present in Fig. ||, the difference Cp^m — Cv,m between the molar 
constant-pressure specific heat Cp^m and the molar constant- volume specific heat Cv,m for 
various pressures. At high temperature Cp^m — Cv,m increases with increasing pressure, while 
at low temperature, it is the contrary and the curve presents a bump. This is can be deduced 



from the following relation BO 



Cp^m — Cv,m — — , (21) 

where ht = ^/Bt- Similarly, the difference Bs — Bt between Bs the bulk modulus calcu- 
lated from F at constant entropy and Bt the bulk modulus calculated from F at constant 
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temperature presents the same behaviour as Cp^m — Cv,m as shown in Fig.|. Indeed, we can 
write ra: 



i^s-^T = (22) 

where ks = l/Bs- 



III. ATOMIC TEMPERATURE FACTOR 

At finite temperature T, the intensity of X-ray diffraction from the crystal is reduced, 
due to atomic motion. The atomic temperature factor e~^^'^^ characterizes the oscillations 



of atom K around its equilibrium position. It is defined |31| by: 



e-^^'^) = exp l^-i 5a/3(/^)G„G^ j (23) 

where Ga is the component of scattering wavevector G which is a reciprocal lattice vector, 
and Bap{hi) is the mean square atomic displacement of atom n along the directions a 
and (3: 

= E ^oth ^^^e.(«:|q, /)e*(/.|q, /) (24) 

where is the mass of atom n, and ej(K|q, /) is the ith component of the eigenvector 
associated with the mode 1 at q in the lattice coordinates. 

In the case of silicon, there is only one kind of atom. Thus, for all atoms e"^*^*^^ are 
identical. The reduction of the diffusion intensity is given by e"^'^^'') , which is usually called 
the Debye- Waller factor. The symmetry of the diamond structure imposes that: 

B^^{k) = B5^p. (25) 

(This B is not to be confused with the bulk modulus.) 

In Table 0, the results obtained for the mean square atomic displacement, at three 
different volumes, corresponding to lattice constants a of 10.00, 10.18, and 10.26 Bohr re- 
spectively, are compared with experimental results. The agreement is on the order of a few 
percent. 
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Interestingly, even at room temperature, the mean square atomic displacement decreases 
with increasing volume, contrary to the intuition. This can easily be understood from the 
definition of B (see Eq. ^). Indeed, it is written as a sum over all phonon bands of a 
term which is proportional to the inverse square of the frequency of the mode u;"^ (since 
cothx for small x). Thus, it is determined mostly by acoustic branches. This 

is confirmed by a band-by-band decomposition of B (see Table where the two lowest 
bands account for more than 2/3 of -B. As the TA band and the first LA band exhibit 
a negative mode Griineisen parameter, we see that their contribution to B decreases with 
increasing volume, at the contrary of the contribution of the second LA band and optic 
bands. 

The thermal parameter B and the atomic temperature factor e~^^'^\ for diffraction with 
scattering vector G = (27r/ao)z, are also calculated as a function of temperature (see Figs. |0 
and 0). The atomic temperature factor is not 1 even at K due to the zero-point motion. 



CONCLUSIONS 

In this paper, dynamical properties of silicon have been calculated using a variational ap- 
proach to density-functional perturbation theory. We have presented first an ab initio study 
of the volume dependence of interatomic force constants up to twenty- fifth nearest neighbors. 
Phonon frequencies of TA(X) and TA(L) modes, and of the associated mode Griineisen pa- 
rameters have also been calculated for different volumes. The influence of successive nearest 
neighbors shells has been analysed. This study has conflrmed that the contribution of atoms 
connected by zig-zag chain along [1 1 0] direction is the most important. It has also proven 
that the contributions of flfth, sixth, and seventh atoms along the chain (respectively thir- 
teenth, eightenth and twenty-flfth nearest neighbors) are not negligible. Analytical formulas, 
taking into account interactions up to second nearest neighbors, have been developped for 
phonon frequencies of TA(X) and TA(L) modes and the corresponding mode Griineisen 
parameters. The volume and pressure dependence of various thermodynamic properties 
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(specific heat, bulk modulus, thermal expansion) were also analyzed. We have pointed out 
the effect of the negative mode Griineisen parameters of the acoustic branches on these 
properties. The effect of zero-point motion was also investigated. Finally, we have presented 
the evolution of the mean square atomic displacement and of the atomic temperature factor 
with the temperature for different volumes, emphasizing the anomalous effects due to the 
negative mode Griineisen parameters, present at all investigated temperatures. 
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FIGURES 

FIG. 1. Vibrational motion corresponding to the (a) TA{X) mode and (b) TA(L) mode in 
silicon. The displacements of ions are along the [1 1 0] direction in (a) and along the [1 1 2] 
direction in (b). 

FIG. 2. Forces induced (linear response) along the [1 1 0] zig-zag chain by the unit displacement 
of a generic atom (in grey) of this chain in the (a) [1 1 0], (b) [0 1], and (c) [1 1 0] directions. 
The unit vectors corresponding to these directions are respectively x, yjand z. The arrows starting 
from white atoms represent the direction of the forces, while the absolute value of it is written 
close to it. Forces along the [1 1 0] direction are indicated by dots, whereas crosses refer to forces 
in the opposite direction. The forces are expressed in Hartree/Bohr. 

FIG. 3. Volume dependence of the Helmholtz free energy F for 4 different temperatures. The 
smallest is K (upper curve), the biggest is 300 K (lower curve) with an increment of 100 K between 
the curves. The equilibrium volume Vo(T) is located at the minimum of each curve. Between 20 K 
and 120 K, this volume decreases due to the negative thermal-expansion coefficient. 

FIG. 4. Temperature dependence of the theoretical bulk modulus Bt for 4 different pressures. 
The smallest is Pa (lower curve), the biggest is 6 10^ Pa (upper curve) with an increment of 2 
10^ Pa between the curves. Bt is expressed in Mbar. Temperature is in K. 

FIG. 5. Temperature dependence of the constant-pressure specific heat Cp^m for 2 different 
pressures. The smallest is Pa (upper curve at high temperature), the biggest is 6 10^ Pa (lower 
curve at high temperature). Cp^m is expressed in J mol^^ K^^. The crosses indicate experimental 
data [28|. Temperature is in K. 

FIG. 6. Temperature dependence of the volumic thermal-expansion coefficient ap for 4 different 
pressures. The smallest is Pa (upper curve), the biggest is 6 10^ Pa (lower curve) with an 
increment of 2 10^ Pa between the curves. Up is expressed in K~^. The crosses indicate experimental 
data [27|. Temperature is in K. 
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FIG. 7. Temperature dependence of overall Griineisen parameter 7 for 4 different pressures. 
The smallest is Pa (upper curve), the biggest is 6 10^ Pa (lower curve) with an increment of 2 
10^ Pa between the curves. The crosses indicate experimental data [^2|. Temperature is in K. 

FIG. 8. Temperature dependence of the difference between the constant-pressure specific heat 
Cp^rn and the constant-volume specific heat Cv,m for 4 different pressures. The smallest is Pa 
(lower curve at high temperature), the biggest is 6 10^ Pa (upper curve at high temperature) with 
an increment of 2 10^ Pa between the curves. Cp,m is expressed in J mol^^ K^^. 

FIG. 9. Temperature dependence of the difference between Bs the bulk modulus calculated 
from F at constant entropy and Bt the bulk modulus calculated from F at constant temperature, 
for 4 different pressures. The smallest is Pa (lower curve at high temperature), the biggest is 
6 10^ Pa (upper curve at high temperature) with an increment of 2 10^ Pa between the curves. 
Bs — Bt is expressed in Mbar. Temperature is in K. 

FIG. 10. Temperature dependence of mean square atomic displacement B for silicon atoms in 
bulk silicon at three different volumes, corresponding to lattice constants a of 10.00 (dashed line), 
10.18 (solid line), and 10.26 Bohr (dotted line) respectively. The values of B are expressed in A^. 
Temperature is in K. 

FIG. 11. Temperature dependence of atomic temperature factor e~^^'^^ for silicon atoms in 
bulk silicon for diffraction with the scattering vector G = (27r/ao)z at three different volumes, 
corresponding to lattice constants a of 10.00 (dashed line), 10.18 (solid line), and 10.26 Bohr 
(dotted line) respectively. Temperature is in K. 
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TABLES 

TABLE L Interatomic force constant matrix elements and associated force Griineisen parame- 
ters for silicon at three different volumes, corresponding to lattice constants a of 10.00, 10.18, and 
10.26 Bohr respectively. The convention of Ref. for labelling the matrix elements has been 
followed. One atom is at the origin, while the coordinates of the second atom are expressed in 
units of a/4, the stars indicate the atoms belonging to the zig-zag chain along [1 1 0] direction. 
The first column indicates the number of the shell A^A^ to which the second atom belongs. The 
force Griineisen parameters of IFC matrix elements that are smaller than 10^'* Hartree/Bohr^ are 
not mentioned. The interatomic force constant matrix elements are expressed in Hartree/Bohr^. 
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Coordinate 




a=10.00 


a=10.18 


a=10.26 






a=10.00 


a=10.18 


a=10.26 





(0,0,0)* 


ao 


0.15629 


0.13904 


0.13201 


7 


'«o) 
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(1,1,1)* 


ai 
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-0.03385 
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^Q!l) 
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7 
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TABLE II. Phonon frequencies of TA{X) mode and associated mode Griineisen parameters for 
silicon at three different volumes, corresponding to lattice constants a of 10.00, 10.18, and 10.26 
Bohr respectively. The reference value (r.v.) is obtained by taking into account interactions with 
all the atoms included in the real space box defined by our sampling of the Brillouin zone (10 
special points). The other values are obtained by limiting the interactions to the successive nearest 
neighbors (NN) shells. The phonon frequencies are expressed in cm~^. 



NN 


a=10.00 


^TA{X) 

a=10.18 


a=10.26 


a=10.00 


7ta(x) 
a=10.18 


a=10.26 





383.462 


361.688 


352.426 


1.080 


1.102 


1.107 


1 


196.622 


206.021 


209.256 


-1.016 


-0.727 


-0.598 


2 


148.056 


161.684 


166.347 


-1.961 


-1.340 


-1.082 


3 


135.107 


150.376 


155.576 


-2.413 


-1.609 


-1.287 


4 


131.783 


147.576 


152.949 


-2.560 


-1.695 


-1.352 


5 


106.607 


128.275 


135.460 


-4.389 


-2.633 


-2.016 


6 


105.327 


127.507 


134.829 


-4.557 


-2.705 


-2.060 


7 


105.390 


127.438 


134.725 


-4.525 


-2.692 


-2.053 


8 


122.392 


142.025 


148.656 


-3.438 


-2.180 


-1.711 


13 


128.111 


146.863 


153.222 


-3.131 


-2.019 


-1.595 


18 


122.443 


142.249 


148.906 


-3.473 


-2.189 


-1.710 


25 


119.535 


140.618 


147.431 


-3.690 


-2.280 


-1.762 


r.v. 


119.903 


140.466 


147.347 


-3.690 


-2.295 


-1.782 
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TABLE III. Phonon frequencies of TA(L) mode and associated mode Griineisen parameters 
for silicon at three different volumes, corresponding to lattice constants a of 10.00, 10.18, and 10.26 
Bohr respectively. The reference value (r.v.) is obtained by taking into account interactions with 
all the atoms included in the real space box defined by our sampling of the Brillouin zone (10 
special points). The other values are obtained by limiting the interactions to the successive nearest 
neighbors (NN) shells. The phonon frequencies are expressed in cm~^. 



NN 


a=10.00 


a=10.18 


a=10.26 


a=10.00 


7rA(L) 
a=10.18 


a=10.26 





383.462 


361.688 


352.426 


1.080 


1.102 


1.107 


1 


147.432 


151.460 


152.695 


-0.609 


-0.396 


-0.295 


2 


119.157 


127.508 


130.363 


-1.494 


-1.041 


-0.845 


3 


120.881 


129.204 


132.089 


-1.459 


-1.033 


-0.848 


4 


124.496 


132.393 


135.122 


-1.346 


-0.954 


-0.783 


5 


90.610 


103.339 


107.667 


-3.001 


-1.951 


-1.547 


6 


79.247 


93.544 


98.304 


-3.889 


-2.390 


-1.845 


7 


76.339 


91.167 


96.078 


-4.197 


-2.534 


-1.942 


8 


98.491 


110.644 


114.798 


-2.630 


-1.745 


-1.396 


13 


101.679 


113.235 


117.210 


-2.416 


-1.628 


-1.312 


18 


96.809 


108.985 


113.149 


-2.680 


-1.776 


-1.420 


25 


96.631 


108.918 


113.746 


-2.532 


-1.959 


-1.737 


r.v. 


96.239 


108.626 


112.867 


-2.742 


-1.814 


-1.451 
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TABLE IV. Mean square atomic displacement B of silicon atoms at T=295 K for three different 
volumes (lattice constants a of 10.00, 10.18, and 10.26 Bohr), and the corresponding experimental 
data at a=10.26 Bohr. The lower part of the table presents a band-by-band decomposition of B. 
The values are expressed in A^. 





Present work 
a=10.00 a=10.18 a=10.26 


Experimental data (a=10.26) 
Ref. H Ref. |34] Ref. || Ref. 


36| 


B 


0.4967 0.4745 0.4707 


0.4613 0.4500 0.4515 0.4660 


B{TA) 
5(LAi) 
5(LA2) 
-B(LO+TO) 


0.2298 0.2052 0.1992 
0.1592 0.1510 0.1482 
0.0496 0.0543 0.0563 
0.0582 0.0641 0.0672 
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